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A key factor influencing a drug’s efficacy is its residence time in the binding pocket of the host 
protein. Using atomistic computer simulation to predict this residence time and the associated 
dissociation process is a desirable but extremely difficult task due to the long timescales involved. 
This gets further complicated by the presence of biophysical factors such as steric and solvation 
effects. In this work, we perform molecular dynamics (MD) simulations of the unbinding of a 
popular prototypical hydrophobic cavity-ligand system using a metadynamics based approach that 
allows direct assessment of kinetic pathways and parameters. When constrained to move in an 
axial manner, we find the unbinding time to be on the order of 4000 sec. In accordance with 
previous studies, we find that the ligand must pass through a region of sharp dewetting transition 
manifested by sudden and high fluctuations in solvent density in the cavity. When we remove 
the steric constraints on ligand, the unbinding happens predominantly by an alternate pathway, 
where the unbinding becomes 20 times faster, and the sharp dewetting transition instead becomes 
continuous. We validate the unbinding timescales from metadynamics through a Poisson analysis, 
and by comparison through detailed balance to binding timescale estimates from unbiased MD. 
This work demonstrates that enhanced sampling can be used to perform explicit solvent molecular 
dynamics studies at timescales previously unattainable, obtaining direct and reliable pictures of the 
underlying physio-chemical factors including free energies and rate constants. 


I. INTRODUCTION 

The unbinding of ligands from host substrates is a phe¬ 
nomenon widely occurring across biological and chem¬ 
ical sciences. It is of great interest to be able to un¬ 
derstand the thermodynamics and kinetics of such pro¬ 
cesses, especially how they are influenced by solvent and 
steric effects. An accurate estimate of unbinding ki¬ 
netics is in fact of crucial importance for drug discov¬ 
ery paradigms [1] [2 . However, in spite of the advent of 
massively parallel computer resources, it has not been 
so easy to simulate the dynamics of ligand unbinding 
and calculate associated rate constants. The compli¬ 
cations are mainly two fold. First, as has been seen 
in studies of model systems [3HT0]. various proteins [IT 
[13] . human immunodeficiency virus [14] and actual anti¬ 
cancer drugs p~5RT7] . the solvent often manifests itself 
at the molecular scale. While coarse-grained models 
can be fit to explicit solvent molecular dynamics (MD) 
simulations 0, predictive power can be attained only by 
performing all-atom MD. Second, performing all-atom 
MD for unbinding of such systems is however plagued 
by the timescale problem. MD is restricted to integra¬ 
tion timesteps of a few femtoseconds, which can be par¬ 
tially mitigated by multiple timestep MD algorithms [18] . 
However, it is not yet routinely feasible to go into the 
millisecond regime and beyond for any system with more 
than few thousand atoms. 

In this paper we consider a popular prototypical cavity- 
ligand system in explicit water where the attraction be¬ 
tween water and the two nanoscale objects, namely a 
fullerene molecule and a spherical cavity, is weak lamier 


unj. We provide a full dynamical picture of the unbinding 
process demonstrating the clear role of water. We find 
that even in this relatively simple system there exists a 
rich range of dynamics that changes qualitatively and 
quantitatively as a function of the cavity-ligand distance 
and the motional degrees of freedom. Previous pioneer¬ 
ing studies 00 uni H2H21] involving explicit all-atom 
MD, brownian dynamics, transition path sampling and 
other approaches have clearly shown that the association 
in such systems has a clear signature of solvent fluctua¬ 
tions and a sharp dewetting transition as the nanoscale 
objects approach each other - if sterically constrained to 
move along the axis of symmetry. This is a popular set¬ 
up that has been considered in numerous studies over the 
years [3. EHEl 10] . and is suggestive of biological systems 
where steric hindrances in the binding pocket do not al¬ 
low the ligand to roll or move in a free manner [121 [Ill- 
Using unbiased MD and brownian dynamics tools, it was 
previously possible to calculate the timescales of associa¬ 
tion or binding for such systems that explicitly accounted 
for the dewetting transition[3 . 

However, apart from one recent work m, to the best 
of our knowledge there is no reported study in which the 
timescales of the analogous dissociation or unbinding pro¬ 
cess were calculated through MD simulations. For such 
cavity-ligand systems due to the very high energy barrier 
of ^30—40 &bF[ 3], the unbinding timescales are simply 
too slow to be amenable through unbiased MD calcula¬ 
tions. As such, in this work we use the popular enhanced 
sampling technique metadvnamics [22H25] along with its 
recent extension [12, 26] for obtaining unbiased dynamics 
to calculate free energy profiles and unbinding rate con- 
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stants for the cavity-ligand system. Furthermore, we ask 
and answer the following question: how does the dynam¬ 
ics of cavity-ligand association and dissociation depend 
on motional degrees of freedom? That is, how would the 
timescales of binding/unbinding vary between the cases 
where the ligand can/can not undergo free to move mo¬ 
tion in any direction? 

We find that the unbinding, in the case of motion being 
sterically constrained along the axis of symmetry of the 
system, proceeds through a sharp wetting transition at a 
critical ligand-cavity separation, conforming to the pic¬ 
ture presented by Mondal et al through their position de¬ 
pendent friction calculation [3]. The mean unbinding time 
of this system is found to be around 4000 seconds. The 
transition pathways harvested from our metadynamics- 
assisted MD runs are in perfect accordance with the pre¬ 
vious calculations of Mondal et al j3]. However, when the 
motion restraint is removed and the ligand is free to move 
in any direction, it finds an alternate pathway wherein 
there is no abrupt wetting transition, and the mean un¬ 
binding time reduces twenty-fold to 200 seconds. The 
binding times are also reduced. We validate rigorously all 
free energy profiles and rate constants through a combi¬ 
nation of umbrella sampling, unbiased MD when feasible, 
and principle of detailed balance. The rate constants are 
further validated also through the Kolmogorov-Smirnov 
test for Poisson distribution proposed in Ref. [27], thus 
giving further confidence in the estimated dynamics. 

This work thus provides useful new insight into 
the phenomenon of hydrophobic interaction in solvated 
nanoscale systems [28], allowing one to directly simulate 
the unbinding process in MD in spite of the very high as¬ 
sociated barriers. It also demonstrates that with a careful 
use of recently developed enhanced sampling techniques, 
one can perform molecular dynamics studies of unbind¬ 
ing/binding that extend well into the seconds timescale, 
and provide statistically accurate thermodynamic and ki¬ 
netic information. The ideas that are used in this work 
are fairly generic and should be applicable to a large 
range of studies pertaining especially to ligand unbinding 
in explicit solvent. 


ii. METHODS 

Throughout this work, we consider a cavity-ligand sys¬ 
tem (Fig. [l]) where the attraction between water and 
both the nanoscale objects is weak and leaning towards 
a hydrophobic system (see SI for detailed potential forms 
and parameters). While the primary method in this pa¬ 
per is metadynamics, we directly and indirectly validate 
any findings from metadynamics with alternate indepen¬ 
dent approaches. 


A. Metadynamics for free energy reconstruction 

Metadynamics is a widely used method for exploring 
complex free energy surfaces characterized by high free 
energy barriers [22lT25| [29] . One first identifies a small 
number of slowly changing order parameters, called col¬ 
lective variables (CVs) [50]. A memory dependent bi¬ 
asing potential is constructed through the simulation as 
a function of these CVs, typically in the form of repul¬ 
sive Gaussians added wherever the system visits in the 
CV space. Thus the system slowly starts to avoid the 
places where it has already visited. This leads to a grad¬ 
ual enhancement of the fluctuations in the CVs, through 
which the system is discouraged from getting trapped 
in the low free energy basins. At the end of a metady¬ 
namics run the probability distribution of any observable, 
whether biased directly or not, can be computed through 
a reweighting procedure ]24l I5T ] . This easy reweighting 
functionality is one of the many features of metadynam¬ 
ics that has made it a very popular method for calculation 
of free energy surfaces. 



FIG. 1: Cavity-ligand system in explicit water with axes 
marked. Red: fullerene shaped ligand atoms. Orange: 
cavity atoms that interact with the ligand. Green: wall 
atoms. See SI for corresponding interaction potentials. 

For the free energy reconstruction through metady¬ 
namics, in the case when the ligand is sterically con¬ 
strained to move along the axis of symmetry (z-axis in 
Fig. 0 , we perform one-dimensional metadynamics with 
the ^-coordinate as the only CV. For the case when the 
ligand is free to move in any direction, we perform two- 
dimensional metadynamics with z- and the radial dis¬ 
tance (p = yjx 1 + y 2 ) from the axis of symmetry as the 
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two CVs. In either case bias is added every 600 fem¬ 
toseconds. In the SI we report the values of all other 
relevant parameters for both cases. In either case, we 
use restraining walls at high ligand-cavity separation to 
facilitate multiple re-entry events (see SI for details of 
restraints). 

B. From metadynamics to dynamics 

Recently, Tiwary and co-workers extended the scope 
of metadynamics by showing how to extract unbiased 
rates from biased ones with minimal extra computational 
burden|26]. For this, they made two key assumptions on 
the dynamics: 

1. the transition processes are characterized by move¬ 
ments from one stable state to another via dynam¬ 
ical bottlenecks that are crossed only rarely, but 
when such a transition does happen, the time spent 
in the bottleneck is small 

2. while there is no need to know beforehand the pre¬ 
cise nature or location of these bottlenecks, one 
should have CVs that can distinguish between all 
stable basins of relevance. Note that this CV does 
not have to be the true reaction co-ordinate [26l 132] . 

Under these two key assumptions, by making the bias 
deposition slower than the time spent in dynamical bot¬ 
tlenecks, it is possible to keep the transition states (TS) 
relatively bias-free through the course of metadynam¬ 
ics. This so-called infrequent metadynamics approach 
m [261 [27] preserves the unbiased sequence of state to 
state transitions and allows one to access the acceleration 
of transition rates achieved through biasing by appealing 
to generalized transition state theory [33] and calculating 
the following simple running average [261, l34l [35] : 

a = (e 0v ^’% ( 1 ) 

where s is the collective variable being biased, (3 is the in¬ 
verse temperature, V(s, t ) is the bias experienced at time 
t and the subscript t indicates averaging under the time- 
dependent potential. The above expression is valid even 
if there are multiple intermediate states and numerous 
alternative reactive pathways [26l [27] . 

In a successive work, a way was also proposed to assess 
the reliability of the two assumptions above[27], This 
relies on the fact that the escape times from a long- 
lived metastable state obey a time-homogeneous Poisson 
statistics[27] with a single rate law. A statistical analysis 
based on the Kolmogorov-Smirnov (KS) test can quanti¬ 
tatively assess how precisely the above assumptions have 
been met m- Thus, if (a) significant bias got deposited 
in the TS region even with infrequent biasing, or (b) 
there are hidden unidentified timescales at play that the 
CV does not resolve, it would lead to failing the test for 
time-homogeneous Poisson statistics. 


For the estimation of kinetics through metadynam¬ 
ics for the sterically constrained case, we perform one¬ 
dimensional metadynamics with the ^-coordinate as the 
only CV (Fig. [I]), but with a relatively much slower bias 
deposition rate of once every 10 ps. For the case when 
the ligand is free to move in any direction, we perform 
one-dimensional metadynamics with the overall ligand- 
cavity separation (i.e. d = y/x 2 + y 2 + z 2 ) as CV and 
bias added every 10 ps. For each case, 14 independent 
simulations were started with the ligand fully bound to 
the cavity. The simulations were stopped when the ligand 
was fully unbound for the first time (see Results section 
for precise definition of unbound), and the respective un¬ 
binding times were calculated using the acceleration fac¬ 
tor (Eq. [l]). The transition time statistics so obtained 
was then subjected to a Poisson analysis to ascertain its 
reliability [27]. In the SI we report the values of all other 
relevant parameters for various cases. 

One of the main features of infrequent metadynamics 
is that the segments of trajectories that cross the barrier 
between successive bias depositions are representatives 
of the unbiased transition path ensemble [36]. As such 
we also present typical reactive trajectories for the sharp 
dewetting transition when it happens. 

As a further verification of the unbinding timescale, we 
back-calculate the corresponding binding timescale using 
the free energy difference between the bound and un¬ 
bound states and the principle of detailed balance. We 
compare this with the estimate of binding time from sep¬ 
arate unbiased MD runs, given that binding is tractable 
through the latter. Specific details of this protocol are 
provided in the respective Results sections. 


C. Umbrella sampling and related fixed bias 
method 

To further examine the quality of the free energies 
obtained from metadynamics, we perform independent 
umbrella-sampling based free energy calculations. For 
the scenario of ligand sterically constrained to move along 
the axis of symmetry, we use the method previously de¬ 
scribed by Mondal et al. [3]. Briefly, in this method, one 
restrains the position of ligand along z direction (Fig. 
|T]) at a given ligand-pocket separation and computes the 
mean force being experienced by the ligand. By doing a 
scan of mean force calculation along a wide-range of val¬ 
ues of z (z— 0.8 to 2.0 nm), we obtain a profile of mean 
force acting between ligand and cavity along the z di¬ 
rection. Finally, we integrate this mean-force-profile to 
obtain the free energy profile along ligand-pocket sepa¬ 
ration z. The solvent-induced part of this free energy 
profile obtained in this fashion was already reported in 
our previous work[3] and here we use the total free energy 
profile for comparison with metadynamics estimate. 

For the scenario where there is no constraint on the 
movement of the ligand, we map the free energy using 
two-dimensional umbrella sampling. For this, we employ 
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FIG. 2: Free energy profiles when ligand is sterically constrained to move along system’s axis of symmetry, (a) 1-d 
free energy as function of z. Blue denotes values from metadynamics with error bars and red from umbrella 
sampling, (b) 2-d free energy as function of 2 and pocket water occupancy w (see [3] for precise definition). The 2-d 
free energy was obtained by re weighting [24] 1-d metadynamics run performed along z. The solid black vertical line 
marks the critical z c value from Ref. [3] at which dewetting transition had been predicted using a completely 
different approach. All energies are in kJ/mol. z is in nm. 

the same two CVs, namely the z- and p = yjx 2 + y 2 
components of cavity-ligand distance, which we use in 
the corresponding metadynamics part as well (subsection 
“Metadynamics for free energy reconstruction”). See SI 
and [3] for further details of the umbrella sampling pro¬ 
tocol. 


hi. RESULTS 

We now systematically present the various results of 
our study. In Table 1, we summarize collectively the 
various timescales for binding and unbinding. 


A. Ligand sterically constrained to move along axis 
of symmetry 

1. Thermodynamic profile 

In Fig. 2(a-b), we report the 1-dimensional free energy 
as a function of the distance £ between the centres of mass 
of the ligand and the substrate. The free energy profile 
as obtained by metadynamics is in very good agreement 
with that from umbrella sampling. We find an energy 
barrier in unbinding of around 80 kJ/mol. At an interme¬ 
diate pocket-ligand separation, there is a relatively small 
but non-negligible barrier to binding as well, of around 
20 kJ/mol. Our close inspection of the successful reactive 
trajectories (Fig. 3) shows that at this location there are 
large fluctuations in cavity-water density that must hap¬ 
pen before the ligand unbinds from the host. Our previ¬ 
ous prediction of dewetting-transition-mediated cavity- 
ligand recognition is confirmed by the two-dimensional 


free energy surface as a function of the ligand-host sepa¬ 
ration z and the number of water molecules in the binding 
pocket, as depicted in Fig. 2(b). Fig. 2(b) demonstrates 
that the sharp solvent fluctuation happens in a particu¬ 
lar ligand-cavity separation of 1.15-1.25 nm, which is in 
very good agreement with previous estimate by Mondal 
et al |3 using an independent approach. 


2. Kinetics 

Having established that 1-dimensional metadynamics 
performed using z as a CV perfectly reproduces the 
salient thermodynamic features of this system including 
the dewetting transition at the correct ligand-substrate 
separation, we then calculate the timescale of unbinding, 
which is the crux of the current article. For this, we use 
the infrequent metadynamics set-up, and calculate the 
acceleration factor achieved in metadynamics. 14 inde¬ 
pendent simulations were started with the ligand docked 
to the host, and stopped when z = 1.4 was attained. We 
find a mean unbinding time of 3863 db 1032 sec (Table 
1). The transition time statistics so obtained fit very well 
a Poisson distribution [271 (see SI for detailed analysis), 
demonstrating that (1) the bias deposition was infrequent 
enough to not gradually corrupt the transition states, (2) 
biasing the CV z was sufficient to ensure the time-scale 
separation needed for the infrequent metadynamics ap¬ 
proach to be applicable. 

To further validate the unbinding time estimate of 
3863±1032 sec, we invoke the principle of detailed bal¬ 
ance. We use the 1-d free energy profile (Fig. 2) and the 
relation unbinding time = binding time xe“^ AG , where 
AG is the free energy difference between the solvated 
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Unbinding (sec) 


Binding (picosec) 


Sterically constrained free to move Sterically constrained free to move 

Metadynamics: 3863 ± 1032 200 ± 51 769 ± 198 118 ± 31 

Unbiased MD: N.A. N.A. 476 ± 33 157 ± 5 


TABLE I: Binding and unbinding timescales in seconds and picoseconds respectively. The binding timescales from 
metadynamics denote values indirectly obtained after use of detailed balance with AG estimate from metadynamics. 



FIG. 3: A typical reactive trajectory obtained via 
infrequent metadynamics providing a clear description 
of the dewetting transition as the sterically constrained 
ligand crosses over from bound to unbound state. Each 
snapshot proceeding row-by-row from top left to 
bottom right represents 10 ps long MD trajectory 
projected on (z,w) space, x-axis here is the ^-distance, 
and y-axis is re, the number of pocket waters. The solid 
black vertical line marks the critical z -value at which 
Ref. [3] had predicted the existence of large solvent 
fluctuations. The underlying free energy surface, axes 
and the contours are same as in Fig. 2(b). 


and bound ligand. Here the mean binding time means 
the time for the fully solvated ligand (z = 1.4) to get 
bound. Note that for the purpose of checking detailed 
balance, as is our objective here, any 2 :—value would be 
fine. As per Fig. 2, the estimate of A G= -72.5 kJ/mol 
from metadynamics, giving us a mean unbinding time of 
769 zb 198 ps, which is well within order of magnitude 
agreement with the value of 476 ± 33 ps through unbi¬ 
ased MD simulations reported by Mondal et al [3]. 

In Fig. 3 we provide a set of snapshots of a typi¬ 


cal ^60ps long reactive trajectory as the system moves 
successfully from bound to unbound state. Here each 
snapshot corresponds to 10 ps of MD projected on the 2- 
dimensional (z,w) space, where w is the number of pocket 
waters (see [3] for precise definition of w). This pro¬ 
vides extremely clear dynamical evidence of the dewet¬ 
ting transition in such a hydrophobic and sterically con¬ 
strained system. The solid black vertical line in Fig. 3 
marks the critical z-value at which Ref. 3 had predicted 
the existence of large solvent fluctuations. Our MD thus 
qualitatively and quantitatively validates the prediction 
of that and previous works. In Table 1, we provide a sum¬ 
mary of the timescales for all the cases obtained through 
metadynamics, unbiased MD and through the use of de¬ 
tailed balance. 


B. Ligand free to move in any direction 

1. Thermodynamic profile 


Fig. 4(a) shows the 2-d free energy as a function of z 
and p = y/x 2 + y 2 obtained through metadynamics. For 
all practical purposes this free energy surface is indistin¬ 
guishable from the one obtained through umbrella sam¬ 
pling and reported in SI. In SI we also provide a compari¬ 
son of the 1-d free energy as a function of 2 obtained from 
metadynamics and from umbrella sampling. Note from 
Fig. 4 that when constraints are lifted, the free energy 
minimum is now slightly off-center and deeper than the 
free energy at the geometric center of the cavity. We find 
that for this case the system avoids the central dewetting 
pathway and instead takes an alternate route, rolling in 
and out along the sides of the cavity. This pathway is fa¬ 
vorable when compared with the axial symmetric path as 
the number of hydrophobic contacts between the ligand 
and the cavity are optimized. Indeed, even though the 
minimum is now deeper than the sterically constrained 
case, the overall barrier is smaller (Fig. 4 and Fig. 1(b) 
in SI), and as we show in next subsection Kinetics, the 
unbinding is faster. 

Fig. 4(b) shows the 2-d free energy but now as a func¬ 
tion of ligand-cavity separation d = \Jx 2 + y 2 + z 2 and 
the pocket water occupancy w. This is to be contrasted 
with Fig. 2 for the sterically constrained motion case. 
This clearly demonstrates that the water occupancy in 
this case changes gradually without any sharp transitions 
at specific cavity-ligand separations. 
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FIG. 4: Various free energy profiles from metadynamics 
for the case when the ligand is free to move in any 
direction, (a) is as function of (z,p), while (b) is as a 
function of (d,w). Note the clear lack of a sharp 
dewetting transition in (b), in contrast with Figs. 2 and 
3. All energies are in kJ/mol and contours are 
separated by 10 kJ/mol. 


2. Kinetics 

To obtain unbinding kinetics for this case, we per¬ 
form 14 independent infrequent 1-d metadynamics bias¬ 
ing d = \/x 1 + y 2 ± z 2 , starting from fully bound ligand 
and stopping when unbound (see Methods and SI for de¬ 
tailed simulation parameters). For the current purpose, 
we take the unbound ligand as 2 = 1.5, p = 0.6 — 0.8. 
The unbinding time now reduces to 200 ±51 sec, with 
again an excellent Poisson fit demonstrating reliability of 
the choice of CV and the timescales so generated. Using 
the free energy difference of -69.8 kJ (Fig 1(b) in SI) be¬ 
tween the bound and unbound states and the principle 
of detailed balance, we obtain a mean binding time of 
118 ± 31 ps. This is again in good agreement with the 
estimate of 157 ± 5 ps through unbiased MD. 


iv. DISCUSSION 

In this paper, we have provided detailed thermody¬ 
namic and kinetic insights into the binding and unbind¬ 
ing of a popular prototypical ligand-substrate system 
through fully atomistic metadynamics-assisted molecular 
dynamics simulations performed in explicit water. Such 
systems are very relevant to the understanding of a range 
of processes across biology, chemistry, and biochemistry. 
We find that in accordance with previous Brownian dy¬ 
namics based studies mmm, the binding of such a sys¬ 
tem proceeds through a marked dewetting transition, but 
only when the system is sterically constrained to move 
along its axis of symmetry. We calculated the unbinding 
time to be around 4000 secs, with an associated bar¬ 
rier of roughly 80 kJ/mol. When the steric constraint 
is removed, the dewetting transition becomes continu¬ 
ous, and the unbinding proceeds through an alternate 
preferred pathway in which the residence time of the lig¬ 
and decreases 20-fold to 200 secs. These extremely long 
unbinding timescales were obtained through the use of 
metadynamics with its recent extension[26], and were val¬ 
idated through alternate simulation techniques and the 
principle of detailed balance. 

Even though the systems considered here are far sim¬ 
pler than an actual complex drug-protein system, we be¬ 
lieve this work has multiple useful implications. Firstly, 
this is one of the first times that such a quantitatively 
insightful study has been carried on realistic ligand- 
substrate systems with very slow unbinding kinetics, and 
validated through a range of simulation techniques. The 
relative simplicity of the system allowed us to validate the 
timescales through detailed balance by performing thor¬ 
ough unbiased MD simulations of the binding, thereby 
giving confidence in the use of metadynamics type tech¬ 
nique for getting unbinding kinetics in more complex 
protein-ligand systems fl2]. Secondly, our work shows 
how the presence of simple steric constraints can heavily 
influence the role played by molecular solvent. We hope 
that our work will serve as a useful addition to the under¬ 
standing of hydrophobic interactions in solvated systems 
in biology and chemistry. 
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